// simply supported beam by Beam2D element
/*
*                ^ force or moment
*        e1      |      e2
*   .____________.______________.
*  /\                           |
*   p1           p2             p3
*/
#include <gtest/gtest.h>
#include <gmock/gmock.h>

#define private public
#define protected public
#include <element.h>
#undef private
#undef protected

using namespace std;

inline void printStdArray6d(const StdArray6d &arr)
{
    for (int i = 0; i < 6; i++)
    {
        cout << arr[i] << ", ";
    }
    cout << endl;
}

inline void printStdArray6d(const StdArray6d* arr)
{
    for (int i = 0; i < 6; i++)
    {
        cout << (*arr)[i] << ", ";
    }
    cout << endl;
}


int main()
{
    // create particles
    Particle* p1 = new Particle(1, 0, 0);
    Particle* p2 = new Particle(2, 5, 0);
    Particle* p3 = new Particle(3, 10, 0);

    // create section
    Rectangle* sect = new Rectangle(1, 2, 4);

    // create material
    LinearElastic* mat = new LinearElastic(1);
    mat->setDensity(10);
    mat->setE(10000);
    mat->setNu(0.31);

    // create elements
    Beam3D* e1 = new Beam3D(1, p1, p2, mat, sect);
    Beam3D* e2 = new Beam3D(1, p2, p3, mat, sect);

    // create DOF
    DOF d1 = {.key={true, true, true, false, false, false},
              .val={0, 0, 0, 0, 0, 0}};
    DOF d3 = {.key={false, true, true, false, false, false},
              .val={0, 0, 0, 0, 0, 0}};
    p1->constraintDof(d1);
    p3->constraintDof(d3);

    // force
    StdArray6d force = {0, 0, 0, 0, 0, 0};
    force[1] = 10;
    // force = {10, 10, 10, 10, 10, 10};

    // solving parameters
    cout << "critical time step: " << e1->calcCriticalTimeStep() << endl;
    double h = 0.01;
    double zeta = 0.5;
    int nStep = ceil(100 / h);

    for (int i = 0; i < nStep; i++)
    {
        if (i == 0)
        {
            p2->setForce(force);
            p1->solveFirstStep(h, zeta);
            p2->solveFirstStep(h, zeta);
            p3->solveFirstStep(h, zeta);
            e1->setParticleForce();
            e2->setParticleForce();
        }
        else
        {
            p1->solve(h, zeta);
            p2->solve(h, zeta);
            p3->solve(h, zeta);
            p1->clearForce();
            p2->clearForce();
            p3->clearForce();
            p2->setForce(force);
            e1->setParticleForce();
            e2->setParticleForce();
        }
        p1->constraintDof(d1);
        p3->constraintDof(d3);
    }
    // results
    cout << "------- particle force ------" << endl;
    cout << "p1: ";
    printStdArray6d(p1->force());
    cout << "p2: ";
    printStdArray6d(p2->force());
    cout << "p3: ";
    printStdArray6d(p3->force());

    cout << "------- particle displace ------" << endl;
    cout << "p1: ";
    printStdArray6d(p1->displace());
    cout << "p2: ";
    printStdArray6d(p2->displace());
    cout << "p3: ";
    printStdArray6d(p3->displace());

    cout << "------- element force ------" << endl;
    cout << "e1, fx, sfy, sfz, mx, my, mz: " << endl;
    cout << "emfi: " << e1->emfi_.fx << ", " << e1->emfi_.sfy << ", ";
    cout << e1->emfi_.sfz << ", " << e1->emfi_.mx << ", " << e1->emfi_.my;
    cout << ", " << e1->emfi_.mz << endl;

    cout << "emfj: " << e1->emfj_.fx << ", " << e1->emfj_.sfy << ", ";
    cout << e1->emfj_.sfz << ", " << e1->emfj_.mx << ", " << e1->emfj_.my;
    cout << ", " << e1->emfj_.mz << endl;

    cout << "\ne2, fx, sfy, sfz, mx, my, mz: " << endl;
    cout << "emfi: " << e2->emfi_.fx << ", " << e2->emfi_.sfy << ", ";
    cout << e2->emfi_.sfz << ", " << e2->emfi_.mx << ", " << e2->emfi_.my;
    cout << ", " << e2->emfi_.mz << endl;

    cout << "emfj: " << e2->emfj_.fx << ", " << e2->emfj_.sfy << ", ";
    cout << e2->emfj_.sfz << ", " << e2->emfj_.mx << ", " << e2->emfj_.my;
    cout << ", " << e2->emfj_.mz << endl;

    delete p1, p2, p3;
    delete mat, sect;
    delete e1, e2;
}

/* result summery
* force = {10, 0, 0, 0, 0, 0};
* e1: emfi = {-10, 0, 0, 0, 0, 0}, emfj = {10, 0, 0, 0, 0, 0}
* e2: emfi = {-1.95863e-09, 0, 0, 0, 0, 0}, emfj = {1.95863e-09, 0, 0, 0, 0, 0}
*
* force = {0, 10, 0, 0, 0, 0};
* e1: emfi = {1.4659e-11, -5, 0, 0, 0, 3.60568e-09},
*     emfj = {-1.4659e-11, 5, 0, 0, 0, -25}
* e2: emfi = {-1.04958e-11, 5, 0, 0, 0, 25},
*     emfj = {1.04958e-11, -5, 0, 0, 0, -1.67897e-09}
*
* force = {0, 0, 10, 0, 0, 0};
* e1: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
* e2: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
*
* force = {0, 0, 0, 10, 0, 0};
* e1: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
* e2: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
*
* force = {0, 0, 0, 0, 10, 0};
* e1: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
* e2: emfi = {0, 0, 0, 0, 0, 0},
*     emfj = {0, 0, 0, 0, 0, 0}
*
* force = {0, 0, 0, 0, 0, 10};
* e1: emfi = {-0, 1, 0, 0, 0, -5.9297e-10},
*     emfj = {0, -1, 0, 0, 0, 5}
* e2: emfi = {0, 1, 0, 0, 0, 5},
*     emfj = {0, -1, 0, 0, 0, -5.9297e-10}
*
* force = {10, 10, 10, 10, 10, 10};
* e1: emfi = {-10, -4.0017, 0, 0, 0, 2.91475e-09},
*     emfj = {10, 4.0017, 0, 0, 0, -19.9915}
* e2: emfi = {-2.14622e-09, 5.9983, 0, 0, 0, 29.9915},
*     emfj = {2.14622e-09, -5.9983, 0, 0, 0, -3.15634e-09}
*/
